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The dynamic separation into phases of high and low baryon density in a heavy ion collision can 
enhance fluctuations of the net rapidity density of baryons compared to model expectations. We 
demonstrate that event-by-event proton and antiproton measurements can be used to observe this 
phenomenon. We then perform real-time lattice simulations to show how these fluctuations arise 
and how they can survive through freeze out. 
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If the QCD phase transition is first order, matter at 
the appropriate temperatures and densities can form a 
mixed phase consisting of plasma droplets in equilibrium 
with a surrounding hadronic fluid. If formed in ion colli- 
sions, this mixed phase can produce large event-by-event 
fluctuations as the system hadronizes M. In particular, 
extraordinary baryon number fluctuations Q can accom- 
pany the first order transition expected at high baryon 
density 0. 

In this paper we explore the dynamics of phase sepa- 
ration in nuclear collisions. The aims of this paper are 
twofold. First, we study the role of baryon number fluc- 
tuations as a probe of the order of the QCD transition. 
We focus on the high baryon density regime, where the- 
ory H and lattice simulations j| suggests that the QCD 
phase transition is first order in a strict thermodynamic 
sense with baryon density as an order parameter. Our 
work may also apply to RHIC collisions, if the low baryon 
density systems produced at the highest energies approx- 
imate a first order transition HQ . Second, we generalize 
techniques from condensed matter physics R to confront 
phase separation in the highly-noncquilibrium context of 
nuclear collisions. Our framework can be used to system- 
atically address other probes as experimental information 
and theoretical understanding evolve. 

To begin, we describe the character of mixed-phase 
baryon fluctuations and show how they can be measured. 
Measurement is not completely straightforward as, e.g., 
neutrons are not easily observed on an event-by-event 
basis. We then formulate a dissipative-hydrodynamic 
model of phase separation and perform numerical sim- 
ulations for that model. 

QCD with two massless flavors can exhibit a first or- 
der transition whose coexistence curve culminates in a 
tricritical point at temperature T c and baryon chemical 
potential fi c For T > T c and p < p c , a second or- 
der phase transition breaks/restores chiral symmetry. If 
the quark masses are sufficiently large, the second order 
transition is replaced by a smooth transformation (since 
chiral symmetry is explicitly broken) . The first order line 
remains, however, with the tricritical point replaced by a 
critical point in the same universality class as a liquid-gas 
transition. 

At RHIC, baryon density may also serve as an approx- 



imate order parameter for the nearly first order tran- 
sition at small net baryon density. Lattice simulations 
HH and general arguments |||],U snow that the baryon 
susceptibility x at /i = can increase suddenly as tem- 
perature is increased near Th ~ 160 MeV, where the chi- 
ral order parameter and, e.g., the energy density change 
sharply. Jumps in the susceptibility commonly accom- 
pany first order transitions. For a liquid-gas transition, 
X = dp /dp, is proportional to the compressibility: steam 
is much more compressible than water. 

Large fluctuations in baryon number occur during 
phase separation in a first order transition. Figure lb 
shows the phase diagram in the T — p plane |J , where p 
is the baryon density. A uniform system quenched into 
the outer parabolic region will separate into droplets at 
the high baryon density p q surrounded by matter at den- 
sity ph- The net baryon number Nb in a sub- volume of 
the system varies depending on the number of droplets 
in the sub-volume. The variance of the baryon number 
Vb = (N B ) — (Nb) 2 can exceed the equilibrium expecta- 
tion by an amount 



AVB~/(l-/)(AiV B ) 2 , 



(1) 



where / is the fraction of the high density phase in the 
sub-volume V and ANb — (p q — Ph)V. In contrast, 
an equilibrium system follows Poisson statistics, so that 
V B = V + V = {N + N), where N, V and N, V are the 
numbers and variances of baryons and antibaryons and 
N B =N-N. 

Experimenters can search for a "super-poissonian" 
variance such as ([!]) by measuring 



£lp — 



V p ^p - (Np + N ¥ ) 
(N p + N ¥ ) 2 



(2) 



where N p and Np 



are the numbers of protons and an- 
tiprotons in a rapidity interval and V p - p is the variance 
of the net proton number N p —Np. This quantity vanishes 
in equilibrium and is related to the more familiar scaled 
variance lu p = (N p + Np)(l + fl p ). Most importantly, f2 p 
is ideal for our application because of the property 



n p = n B ^ VB - {N + N) 



(N + N) 2 



(3) 
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where N and N are the numbers of baryons and an- 
tibaryons - including unseen neutrons and antineutrons 
(the proof follows). The conditions for which (|^) holds 
are met by a range of thermal and Glauber models that 
respect isospin symmetry. Isospin fluctuations can alter 
(|J) near the tricritical point or in the presence of a disori- 
ented chiral condensate, but those effects will be evident 
from pion measurements. 

We demonstrate (|3|) by writing the joint probability 
for N p and N ¥ as J2 N;N p{Np\N)p(Np\N)P(N, N). The 

distribution P(N,N), which determines Qb, is modified 
by phase separation; we make no assumptions about 
its form. We assume that the conditional probabil- 
ity p(N p \N) for measuring N p given N baryons is bi- 
nomial, with q the chance that any individual baryon 
is a proton (see [jHj for notation). We further take 
p(Np\N) for antiprotons to be binomial with the same q. 
These assumptions hold for most thermal and multiple- 
scattering models. The average of the joint distribution 
is (Np + Np) = J2n~nP(^> -^X^ + i 7 ); where the binomial 
averages v = J2n p p( n p\ n ) n p and v = J2n ¥ p( N p\ n ) N p 
yield (N p + Np) =q(N + N). The quantity (N + N) de- 
pends only on P(N,N). Similarly, we find ((N p -N ¥ ) 2 ) = 
q 2 ((N-N) 2 } + q(l-q)(N + N). We combine these mo- 
ments to obtain (^). 

The antiproton contribution to (0, |J) is large only at 
RHIC, where N ¥ /N p ~ 0.6 at ^ = 130A-GeV jjll). 
At the top SPS energy, we estimate p contributions to 
(|f) to be at the few percent level in Au+Au at ^fs = 
17.5 A-GeV, since N ¥ /N p ~ 6% @. The highest baryon 
density - and the greatest potential for observing a first 
order transition - is perhaps at lower energies. 

We remark that Jeon and Koch and Asakawa et al. 
have proposed that hadronization may change the char- 
acter of charge and baryon number fluctuations even in 
the absence of a phase transition jl3| . This effect is essen- 
tially poissonian, however, so it is not clear that it would 
cause fl p to differ from zero, the equilibrium value, or 
that it could be tested without measuring neutrons. The 
effect on charge fluctuations is much more dramatic [fl3| . 

We now turn to describe the process of phase separa- 
tion. To describe the state of the mixed phase, we follow 
the standard condensed matter practice R] and write a 
Ginzburg-Landau free energy / = K(Vp) 2 /2 + /q, where 

/o = -m 2 (p-p c ) 2 /2 + A(p- Pc ) 4 /4 (4) 

describes the excursions of the baryon density p from its 
equilibrium value in the uniform matter. For m 2 cx T c — T 
we find the correct liquid-gas critical exponents. The val- 
ues ph and p q in fig. 1 correspond to the equilibrium den- 
sities at T < T c : ph = p c — Ap and p q = p c + Ap, where 
Ap = \J m 2 /X. The k term describes the droplet surface 
tension. For our / , we compute a = (8nm 3 /9A 2 ) 1 / 2 cx 

To describe the dynamics of the system, we must ac- 
count for the fact that baryon number is conserved. Fur- 



thermore, it is crucial to include dissipation to describe 
this strongly fluctuating system. The simplest equations 
that meet these criteria are: 

dp/dt = M V 2 M , M = /o - *V 2 P ; (5) 

model B in . We illustrate that (|j) describes diffusion 
in a stable liquid by considering fluctuations about the 
equilibrium density p = p/i + <5pfcexp(— ik • x), where 
Spk <C Ph- A system at this density is near the minimum 
of / , so that « fo'{ph)$Pk = 2m 2 5p k . Therefore, 
(||) is standard diffusion equation at linear order in Sp k - 
We identify the baryon diffusion coefficient at ph as D — 
2m 2 M. In general, diffusion drives the system towards 
homogeneity at all density for which /q (p) > 0. 
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FIG. 1. Free energy (a) and phase diagram (b) vs. baryon 
density for (Q). 

Phase separation is most dramatic if the rapid expan- 
sion of the heavy ion system drives the system into the 
unstable region; i.e., the inner parabolic region in fig. lb, 
corresponding to /o'(p) < in fig. la. Droplets form from 
runaway density fluctuations in a process known as spin- 
odal decomposition. We estimate the time scale tr for 
this process by considering the time evolution of small 
fluctuations near p c . We take p = p c + Spk exp(— ik • x) 
to find 

±5p k = m 2 M (k 2 - 4fc 4 ) 5 Pk = ^ (6) 

to linear order in 5 p. Long wavelength disturbances cor- 
responding to < k < \firn 2 /k) grow with time, while 
the surface-tension term stabilizes the shorter wavelength 
modes. The time scale for growth Tk is shortest at 
k R = V(m 2 /2 K ). 

The fastest-growing mode at wavenumber Ur domi- 
nates the early evolution of the system in the unstable 
regime. The time scale for the growth of this mode is 

T R = &e/D (7) 

where £ = n 1 / 2 /m is the correlation length. The values 
of £ and tr determine the time and spatial scales for the 
onset of spinodal decomposition. We plausibly estimate 
the correlation length to be £ ~ 1 fm, roughly the value 
of the inverse sigma mass. For a value of D ~ 8 fm con- 
sistent with calculations in , we find tr ~ 1 fm as [|| . 
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We remark that the large magnitude of D suggested by 
jl7j is consistent with our assumption in (^) that baryon 
diffusion is the dominant transport mode for baryons at 
high density. Our model superficially suggests a slower 
onset of the instability for a substantially smaller value of 
D. However, if D were truly small then it would be nec- 
essary to include transport mechanisms involving convec- 
tion and viscosity. In fact, viscosity must dominate near 
fj, = 0, where diffusion precesses are strictly irrelevant 
& 

To describe nuclear collisions, we extend (^|) to include 
drift due to Bjorken longitudinal flow: 



dp/dr + P /t = MVV, 



(8) 



where r is the proper time and p is given by ^). The 
new drift term forces the average density to decrease as 
(p) oc t , driving the system through the phase coexis- 
tence region. Fluctuations grow when densities approach 
p c (see fig. 1). To derive the drift term, observe that 
(||) follows from baryon current conservation, which more 
gen erally implies <9 M j M — 0. The current is j M = pu 11 + j% 
|]24| , where it M is a fluid velocity that includes a contribu- 
tion from the meson flow, and jd is the diffusion current, 
oc V/x when u = (1,0,0,0). The left and right sides of 
(H) respectively follow from d^pu^) and d^j^ for Bjorken 
flow. 

For times t 3> tr, the system undergoes a nonlinear 
evolution in which droplets merge, reducing their surface 
energy. To study this regime, we write the evolution 
equation (^|) in the dimensionless form: 



9f 
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v 2 (v>- v 3 + v 2 vO 



(9) 



where we use the dimensionless coordinates f = 8t/tr 
and x = x/£. The dimensionless order parameter ip = 
(p — p c )/ Ap equals ±1 when p = p) % q = p c ± Ap. The 
only remaining parameter is e = p c / Ap, which controls 
the strength of the first order transition. Here, we take 
e = 1 corresponding to a strongly first order transition. 
Observe that (^) depends on the temperature and den- 
sity scale only through e. This is an artifact of our very 
simplistic quadratic /o; we will introduce a more realis- 
tic free energy density in later work to study the role of 
temperature in the evolution. 

We solve (||) numerically on a 2+1 dimensional lat- 
tice following Grant et al. pl|. We use a forward Eu- 
ler method to evolve the system in time for a time step 
Af = 0.05. We study the evolution in the transverse 
plane and in the rapidity t)-xt plane, where xt is a carte- 
sian transverse coordinate. The laplacian in the tj-xt 
case is 



V = h Vj_. 

f 2 drf 



(10) 



To treat the higher spatial derivatives we extend the 
next-nearest-neighbor algorithm developed by Oono and 



Puri |16|] and used in [[L5| to account for the asymmetric 
•q-xx lattice. We write 



v 2 v 



1 



2{Axf 



\NN NNNN 



Ip Ip 



(11) 



where the first sum runs over the four nearest neighbors 
(NN) and the second over the four adjacent next-next- 
nearest neighbors (NNNN) . Oono and Puri use the diag- 
onal next-nearest-neighbors instead in ( pd[ ) - a formula- 
tion that requires a symmetric lattice. We take Ax = 1. 
We find that our results are practically indistinguishable 
from NNN results for this spacing on a symmetric 
lattice. To study longitudinal expansion, it suffices to 
replace one coordinate Ax for tA-t]. 




FIG. 2. Order parameter in the transverse plane in the 
absence of expansion. Droplets tend to merge. 

Figures 2 and 3 show 2+1 dimensional numerical sim- 
ulations of (||) in the transverse plane. Only longitudi- 
nal expansion is considered so the coordinates are carte- 
sian with periodic boundary conditions. For comparison, 
fig. 2 shows results in which expansion is neglected by 
omitting the term (e + ip) /f in (|9|) . Expansion shown in 
fig. 3 prevents droplets from merging as in fig. 2. The 
expanding system reaches p c at To = 5 fm. Because this 
is a dissipative system, we must apply thermal noise at 
each lattice site at r to seed phase separation (noise at 
earlier times is dampened). The memory of the initial 
conditions is essentially lost for t — tq > tr. 

We now study the rapidity dependence of baryon num- 
ber fluctuations. Figure 3 shows the computed variance 
for two different initial times and for two rapidity inter- 
vals. The variance is computed from a sample of 5000 
simulated events, each unique due to the thermal noise. 
We see that the super-poissonian fluctuations grow ap- 
preciably by t ~ 2tq . This variance drops as the rapidity 
interval is increased. We find that variance is governed 
by the ratio tq/tr, which compares the expansion and 
droplet-growth time scales. 
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FIG. 3. Order parameter in the transverse plane including 
expansion. Expansion prevents droplets from merging. 

These schematic calculations serve to illustrate the im- 
pact of droplet formation on baryon fluctuations. To ob- 
tain more quantitative predictions, one must use a more 
realistic form of the free energy fg. Our quadratic fo{p) 
strictly applies only near T c and yields compressibilities 
that are equal at ph and p c . This result is unchanged 
if linear and cubic terms are added. In contrast, lattice 
QCD calculations suggest that the compressibility may 
jump across the transition S. The bag model equation 
of state describes a first order transition and predicts a 
jump A(dp/dfi) ~ 2T C 2 for two light flavors, but is not 
analytic in the two-phase region. We will discuss a more 
sophisticated parameterization of fo(p) in future work. 
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FIG. 4. Enhanced variance vs. time for two rapidity win- 
dows, eq (1). 

Turning now to the interpretation of experiments, we 
emphasize that the identification of phase-transition in- 
duced fluctuations requires a systematic comparison of 
fl p in pp, p-nucleus (pA) and nucleus-nucleus (AA) col- 
lisions. Benchmark RQMD and HIJING simulations for 
central Au+Au collisions yield f2 p « in the absence of 
a phase transition in marked contrast with fig. 4. 
However, we have been unable to find experimental re- 
sults on net-proton fluctuations in pp or pA collisions in 
the literature, so it is not clear whether these benchmark 
estimates are reliable. It is therefore important that AA 



experimenters study fluctuations of identified baryons in 
pp and pA collisions. Measurements of fluctuations of 
unidentified charged particles in hadronic collisions 1 19f l 
and strange baryon production in e + e _ collisions 2Q] 
hint of substantial proton-antiproton correlations in pp 
collisions. If it turns out that light and heavy ion fluctu- 
ations are similar, it may be necessary to correlate baryon 
measurements with other signals to extract phase tran- 
sition information, as in pl| . 

Nevertheless, we stress that it is unlikely that super- 
poissonian fluctuations in nucleon-nuclcon (NN) colli- 
sions - if present - result in significant fluctuations in AA 
interactions unless there is a major source of coherence or 
collectivity. If we treat the AA collision as a superposi- 
tion of NN subcollisions, then fl p (AA) = Q P (NN)/N(b), 
where N(b) is the number of participant nucleons. To ob- 
tain a rough upper bound on il p (NN), we take the total 
charge fluctuations measured to be ~ 0.6 in 200 GeV 
pp collisions from Whitmore's review [fl9|| . For Au+Au 
collisions at b < 10 fm, we estimate il p (Au + Au) = 
Qp(NN)/N(0) < 0.01, where we use the wounded nu- 
cleon model to compute N(b) « 59 for b — 10 fm and 
372 for b = 0. RQMD Au+Au simulations for impact 
parameters fall below this bound flT8| , 

We expect Q p to dramatically increase in heavy ion 
systems compared to light ones. In central S+S we expect 
the NN contribution to O p to be below 1%, as implied 
by our wounded nucleon model estimate. Since there is 
no evidence of a phase transition in such light systems 
at AGS or SPS, the appearance of fluctuations at the 
level of fig. 4 in Au+Au would be impressive. But is 
there any source of coherence or collectivity other than a 
phase transition? Gluon junction effects |^2| can lead to 
correlated baryon production in pp, pA and A A collisions. 
This effect is only partially included in RQMD Q. We 
are currently studying how gluon junctions can effect Q p 

B 

In summary, we have studied the phenomenological 
impact of baryon density, a proposed order parameter 
of the putative first order QCD phase transition at high 
baryon density H. We have shown that phase separa- 
tion in the nonequilibrium heavy ion system can lead to 
large baryon fluctuations. These fluctuations are super- 
poissonian and, consequently, can be extracted by mea- 
suring protons alone. For (|J) with (p) oc r _1 , the system 
is unstable only for r < 2.3 To- We extend the calcu- 
lations to much longer times to demonstrate that the 
fluctuations in rapidity survive well past the freeze out 
time, of order 10-30 fm, in accord with 

For sufficiently large To , final state fluctuations can be 
substantial. However, we have seen that a more rapid 
expansion corresponding to smaller To leads to an "infla- 
tion" that prevents the fluctuations from having a large 
impact on the final state. If experiments find that the 
non-poissonian component of fluctuations is small, we 
must use information from flow signals to ascertain the 
degree of this inflation. 

We emphasize that these calculations include diffu- 
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sion, which dampens the fluctuations once the system be- 
comes stable. While diffusion is the primary mechanism 
for dampening fluctuations at high density, viscosity be- 
comes more important at small net baryon density. Sev- 
eral key questions remain: At what energy do heavy ion 
collisions reach a baryon density where the phase tran- 
sition is strongly first order? Is there a residual modifi- 
cation of fluctuations due to the near transition at zero 
baryon density? To what extent does cooling, convection, 
viscosity and collision-geometry alter £l p compared to our 
estimates? Finally, we note that our mixed-phase effect 
may be compensated to some extent by the effect due 
to the difference between fluctuations in a plasma com- 
pared to a hadron gas Nevertheless, the strength of 
the signal in our exploratory calculations invites further 
work. 
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